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We present a new method to calculate Lyapunov exponents of rigid diatomic molecules in three 
dimensions (12iV dimensional phase space). The spectra of Lyapunov exponents are obtained for 
32 rigid diatomic molecules interacting through the Weeks- Chandler- Anderson(WCA) potential for 
various bond length and densities, and compared with those in Y.-H. Shin et al. Our algorithm 
is easy to implement and total CPU time is relatively inexpensive. 



I. INTRODUCTION 



dV 



One of the ways to quantify the dynamical instabil- 
ity of a many-body system is to examine its Lyapunov 
exponents. Lyapunov exponents are a measure of the 
average rate at which nearbytrajectories converge or di- 
verge in the phase space @,H, 0]. A positive Lyapunov 
exponent indicates a divergence between nearby trajec- 
tories, i.e., a high sensitivity to initial conditions. There 
are various methods to compute Lyapunov exponents in 
the literature d, i, 0, S, i] • 

In this paper we propose a new method to calculate 
Lyapunov exponents of the rigid diatomic molecules in 
three dimensions, which does not require periodic rescal- 
ing of the bond length. Rigid diatomic molecules are of- 
ten described by hard dumbbells and there have been sev- 
eral works on the L yap unov instability for the diatomic 
molecular models[l|,M El El El El • 

Our method is based on the algorithms proposed by 
Omelyan et al. [1 61 ] . They proposed the optimized Verlet- 
like algorithms which are derived on the basis of an ex- 
tended decomposition scheme at the presence of a free 
parameter and which are more efficient than the origi- 
nal Vcrlct versions that corresponds to a particular case 
when the introduced parameter is equal to zero. 

In the present paper, we extend Omelyan's algorithms 
to the calculation of Lyapunov exponents, and show that 
our algorithm works efficiently. We compare present re- 
sults with those appeared in Ref.[l], which uses two coor- 
dinates representations to avoid the singularity occurring 
in the equations of motion by combining with the adap- 
tive Runge-Kutta-Fehlberg method of order four. 



II. EQUATIONS OF MOTION 

The position and the orientation of the diatomic 
molecule are specified by a vector to the center of mass, q, 
and a unit vector which points along the molecular axis. 
We denote this unit vector by S. Then, the equations of 
motion are [13] 



S 
L 
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-L x S = uj x S, 



-S x 



dV 



(i) 



where p represents the momentum, M the total mass 
of the molecule, / the moment of inertia, and L the 
angular momentum. As in the Verlet algorithm, we 
split the Hamiltonian for rigid diatomic molecular mo- 
tion into kinetic and potential parts: H = T + V, where 
T = y KI P-P + jjL ■ L and V = V(q, S). The Weeks- 
Chandler- Anderson(WCA) potential with a cutoff length 
r c = 2 1 / 6 is used as the intermolecular potential function. 
It should be noted that the previous work on the sys- 
tem composed of hard dumbells[12] employs more com- 
plicated algorithm because of the discontinuity imposed 
on the interaction potential. In the present work, all the 
numerical results are based on the simulations performed 
on the system composed of 32 diatomic molecules. 

The exact solutions with initial condition xq = 
(qo,Po, So, Lq) are given as follows. 

f . Motion induced by V : 



dV - - - dV 
exp v (t)f = qo, Po - t—\s Q , So, L - tS x —^\s 

dq dS 



2. Motion induced by T : 



exp T (t)f = ( <f + jtPo, Po, S', L 



(3) 



where we use the algorithm in Ref. [l8| to update S 

_ {So + {ujx So) At + ^jg(g ■ So) - |(g ■ J) So]) 
[1 + (wAi/2) 2 ] 



(4) 
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One of great advantages of this algorithm is that we don't 
need to adjust the bond length during the simulations. 
This is also true in Ref.[l|, where the bond length is 
naturally fixed without using any additional constraint. 
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However, the equations of motion in our model are much 
simpler than those used in Ref. [l| , and they are easy to 
implement. 



A split Hamiltonian scheme for the integration, thus, 
is given as 



J 



[g i + 1 ,p i+1 ,S i+1 ,L i + 1 



expv^Ai) exp T (Ai) eocpv(-At) x (q'.p'.S'.L' 



(5) 
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In the actual simulations, however, we use the optimized 
Verlet-like algorithms proposed recently by Omelyan et 
al.llj^to get the Lyapunov exponents. In their paper it 
is shown that the optimized position- Verlet-like(OPV) 
algorithm is a little more efficient than the optimized 
velocity- Verlet-like(OVV) algorithm. Thus, we use the 
OPV algorithm in this paper. The algorithm reads 





= ff(t)-f 
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= K*H 


-vi At/2, 




= &(*) 


+ -/(r»(l-20At, 

TO 


r(t + At) 


= ri(t) 


+ vii At/2, 


v(t + At) 


= VH + 


-/(f(< + A<)KAi, 

TO 



(6) 



where the parameter f » 0.1931833275037836 and 
fj, fn, vi, and v*ii denote the center of positions and their 
velocities of a diatomic molecule. Similarly, we have the 
propagation of the unit vector S and angular momentum 
L from time t to t + At. A time step At = 0.0005 is used 
in our simulations. 



SS = j6(L xS) = 5{uj x S), 

SL = -5(5 x ^4). 

dS 



(7) 



It should be noted that we consider only the first or- 
der, O(At), for updating SS. In our model we have 12 N 
dimensional phase space. In order to compare our re- 
sults with those from lOiV dimensional phase space in 
Ref. [l|, we subtract the additional two degrees of free- 
doms in each diatomic molecule. We have two constraints 
in equations of motion : 



S-SS = 0, 



S-SL + 8S-L = 0. 



(8) 



(9) 



These equations are obtained from two conditions : S ■ S 
= 1 (normalization) and S ■ L = (orthogonality). In 
order to fulfill two constraints SS and SL are replaced 
with, 



SS = SS - (5 • SS)S, 



(10) 



III. LYAPUNOV EXPONENTS 

The initial configuration of the molecules and the sim- 
ulation method in Ref.[l| are used in our simulations. 
The velocities are repeatedly scaled to adjust the required 
temperature from sufficiently high temperature . Once 
the required temperature is obtained, the system is equi- 
librated for 500 time units (10 6 iterations), and data is 
collected for 500 time units(10 6 iterations) to evaluate 
the Lyapunov exponents. Below all quantities are given 
in reduced units. 

In order to get the exponents we use the classical 
method of Benettin at a/.[l9[ refined by Hoover and 
Posch[2(| [U, [22], that requires continuous orthonor- 
malization. The corresponding equations for the varia- 
tion are given as follows. 



SL = SL-{S -SL + SS ■ L)S, 



(11) 
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respectively. Then, 127V dimensional phase space can be 
reduced to lOiV dimensional one. 

We find the characteristic features of the spectra of the 
Lyapunov exponents in our model are the same as those 
in Ref. [l|. In Fig. Q]we compare the positive branch of 
exponents from our calculation with that from Ref. [l[ at 
T=1.0, D=0.3, _B=1.0, where T is temperature, D is the 
density, and B is the bond length, respectively. A; is the 
discrete spectrum of the Lyapunov exponents, and index 
I represents 1, 160, i.e. the half of total number of all 
phase space variables . 

Table U shows the largest Lyapunov exponents Ai, the 
smallest positive Lyapunov exponent A156, and four van- 
ishing exponents (A157, A158, Aisg, and A160) for various 
bond length B at a fixed number density D=0.3, and 
D=0.5, respectively. Temperature is set to T=1.0 for 
both cases. The same trend for the largest Lyapunov 
exponent Ai was already shown in Ref.[l(. 

On the other hand, the dynamics of the tangent vectors 
in the subspace are slightly different. We calculate the 
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FIG. 2: Mean-squared values of the projection at T—1.0, 
D=0.3, B=1.0. 



mean-squared value of the projection of tangent vectors 
5i onto TX subspace which is defined as usual 

(P x , l ) = (V(X)5 l -V(X)S l ). (12) 

FigfJ shows the projection at T=1.0, D=0.'S, 5=1.0, 
where Q denotes x, y, and z components. In our model 
($Xl) ^ or & an( i L P nase space and their tangent space 
SS and SL are not symmetric with respect to the cen- 
ter. A similar asymmetric behavior was also found in 
Ref. [24|. However, we find the Hamiltonian nature of the 
system, i.e., an increase of instability accumulated in one 
subspace is always accompanied with a decrease of insta- 
bility in its conjugate subspace. In Figj2] we also show 
Sq + Lq, which is symmetric with respect to the center, 
for comparison. 

It is interesting to note that the overall patterns are 
symmetric if temperature is low or the bond length is 
small. In Fig|3]we show the projection at T=1.0, i?=0.2, 
and at T=0.1, fl=1.0, respectively. The density is set to 



TABLE I: The largest Lyapunov exponent Ai and the small- 
est positive Lyapunov exponent A 156, and four vanishing ex- 
ponents (A157, A158, A159, and Aieo) for various bond length B 
at a fixed number density D=0.3 (upper) and D=0.5 (lower). 



B 


0.2 


0.4 


0.6 


0.8 


1.0 


Ai 


4.320 


4.574 


4.791 


4.989 


5.179 


Al56 


0.104 


0.170 


0.274 


0.272 


0.229 


A157 


0.025 


0.046 


0.041 


0.034 


0.030 


Al5S 


0.018 


0.014 


0.014 


0.017 


0.013 


A159 


0.010 


0.008 


0.012 


0.012 


0.011 


Aleo 


0.012 


0.009 


0.009 


0.011 


0.007 




B 


0.2 


0.4 


0.6 


0.8 


1.0 


Ai 


6.158 


6.526 


6.557 


6.327 


5.746 


Al56 


0.108 


0.113 


0.052 


0.052 


0.032 


Al57 


0.068 


0.056 


0.028 


0.015 


0.013 


Al58 


0.020 


0.013 


0.015 


0.014 


0.011 


Al59 


0.014 


0.011 


0.010 


0.011 


0.012 


Al60 


0.010 


0.009 


0.007 


0.010 


0.003 
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TABLE II: Comparison of total CPU time(in hours) for the 
simulations on a single node (Intel 3.06GHz CPU). 





This work(A) 


Shin et afs(B) 


A/B 


T= 


1.0, D= 


=0.3, B= 


0.2 


8.8h 


16.5h 


0.53 


T= 


1.0, D- 


=0.3, B= 


1.0 


15.5h 


24.6h 


0.63 


T= 


1.0, D= 


=0.5, B= 


0.2 


17.1h 


28.3h 


0.60 



Z?=0.3 for both cases. It is not clear to us why the pro- 
jection depends on temperature or the bond length, and 
symmetry becomes broken as temperature or the bond 
length is increasing. We speculate that this might come 
from a relation between S and L. Note that the correct 
conjugate momentum for rotation is|17| 

II = IS = L x S, (13) 

not simply L. More detailed analyses will be needed re- 
garding the asymmetric behavior. It can be shown that 
the overall patterns are always symmetric if our coordi- 
nate system is transformed to that of Ref.[l[, i.e., (q, p, 
9, <fi, pe, P4>) (See Fig [4|. Clearly, this shows that our 



coordinate system is consistent with that of Ref . [l[ . 

In Table lU we show comparison of total CPU time for 
our simulations and Shin et aVs on the Lyapunov expo- 
nents. Our code is relatively inexpensive, although we 
need to implement two constraints in Eq.© and Eq.©. 
In Ref.[l[ the coordinate transformation is needed to 
avoid singularity occurring in the equations of motion, 
and a sophisticated integrator is used because the coor- 
dinate transformation cannot be applied to the calcula- 
tion of the Lyapunov exponents. These might consume 
a relatively large CPU time during the simulations. 

In summary, we propose a new approach to calculate 
Lyapunov exponents of the rigid diatomic molecules in 
three dimensions, which does not need a rescaling of the 
bond length, and it is computationally relatively inex- 
pensive. 
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FIG. 4: Same as in Fig. [2] After transforming of the coordi- 
nate system. 
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